Native Plants Analysis

library(vegan)
library(broom)
library(MASS)
library(tidyverse)

## Do Not Run
#wExp <- readRDS("~/RMB/SoilDomestication/Data/wExp.rds")
#wExp %>% 
#  filter(Site == "Arkansas") %>% 
#  ungroup() %>% 
#  saveRDS("~/RMB/SoilDomestication/Data/native_plants_data.rds")
nExp <- readRDS("~/RMB/SoilDomestication/Data/native_plants_data.rds") %>% filter(SampleID != "Arw.22")
tax <- readRDS("~/RMB/SoilDomestication/Data/gg_otus_tax.rds")

long_pcoa <- function(x, samples = "SampleID", otus = "variable", value = "RA", dist = "bray"){
  to_drop <- otus
  metadata <- x %>% 
    dplyr::ungroup() %>% 
    purrr::discard(is.double) %>% 
    dplyr::select_(.dots = paste("-", to_drop)) %>% 
    dplyr::distinct()
  
  wide_table <- x %>% 
    dplyr::select_(samples, otus, value) %>% 
    tidyr::spread_(otus, value, fill = 0)
  
  pc <- vegan::capscale(log2(wide_table[,2:ncol(wide_table)] + 1) ~ 1, dist = dist)
  axes <- dplyr::bind_cols(metadata, dplyr::as_tibble(vegan::scores(pc, choices = c(1:5))$sites))
  eigen_vals <- vegan::eigenvals(pc) / sum(vegan::eigenvals(pc))
  
  return(list(axes = axes, eigen_vals = eigen_vals))
}
raref <- function(x, sampling_depth, value = "value", otus = "variable") {
  rare_values <- data.frame(table(sample(x$`variable`, sampling_depth, replace = T, prob = x$`value`/x$depth)))
  names(rare_values) <- c(otus, "rare_value")
  return(suppressMessages(left_join(x, rare_values, by = otus) %>% replace_na(list(rare_value = 0))))
}
long_adonis <- function(x, samples = "SampleID", otus = "variable", value = "RA", dist = "bray", formula) {
    to_drop <- otus
    metadata <- x %>% 
        dplyr::ungroup() %>% 
        purrr::discard(is.double) %>% 
        dplyr::select_(.dots = paste("-", to_drop)) %>% 
        dplyr::distinct()
  
    wide_table <- x %>% 
        dplyr::select_(samples, otus, value) %>% 
        tidyr::spread_(otus, value, fill = 0)
    permanova <- vegan::adonis(as.formula(paste("wide_table[,2:ncol(wide_table)] ~ ", formula, sep = "")), data = metadata, dist = dist)
    return(permanova)
}

Phyla Stuff

Differential Abundance

safe_glm <- possibly(glm, NA_real_)
safe_glm.nb <- possibly(glm.nb, NA_real_)
nExp2 <- nExp %>% 
  ungroup() %>% 
  mutate(variable = factor(variable)) %>% 
  group_by(SampleID) %>% 
  nest() %>% 
  mutate(rare = map(data, ~raref(., sampling_depth = 15000, value = "value"))) %>% 
  unnest(rare)
pois_glm_plant <- nExp2 %>% 
  filter(Compartment != "Bulk Soil") %>% 
  group_by(variable) %>% 
  filter(sum(rare_value > 0) / n() > 0.1) %>% 
  mutate(host_common_name = relevel(factor(host_common_name), ref = "Rice")) %>% 
  group_by(variable, Compartment) %>% 
  nest() %>%
  mutate(models = map(data, ~suppressWarnings(safe_glm.pois(rare_value ~ host_common_name + offset(log2(depth)), family = poisson, .))))
nb_glm_plant <- nExp2 %>% 
  filter(Compartment != "Bulk Soil") %>% 
  group_by(variable) %>% 
  filter(sum(rare_value > 0) / n() > 0.1) %>% 
  mutate(host_common_name = relevel(factor(host_common_name), ref = "Rice")) %>% 
  group_by(variable, Compartment) %>% 
  nest() %>%
  mutate(models = map(data, ~suppressWarnings(safe_glm.nb(rare_value ~ host_common_name + offset(log2(depth)), .))))
lm_plant <- nExp2 %>% 
  filter(Compartment != "Bulk Soil") %>% 
  group_by(variable) %>% 
  filter(sum(rare_value > 0) / n() > 0.1) %>% 
  mutate(host_common_name = relevel(factor(host_common_name), ref = "Rice")) %>% 
  group_by(variable, Compartment) %>% 
  nest() %>%
  mutate(models = map(data, ~lm(log2((RA * 1000) + 1) ~ host_common_name, .)))